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ABSTRACT 


Interpolation, or resampling coefficients, which are generated 
from low pass filter Fourier Transforms yield more accurate 
resampled values than those obtained using cubic spline tech- 
niques. This is due to the utilization of six data points rather 
than four as currently used In cubic spline analysis. 

After resampling functions are applied to Image data, artifacts 
which are similar to ringing may become pronounced. These 
effects are often present in the original data and the interpola- 
tion merely enhances them. 
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INTERPOLATION ALGORITHMS 
AND IMAGE DATA ARTIFACTS 


I. INTRODUCTION 

During the past few years, a great deal of discussion has taken place at various 
conferences and In various reports regarding the application of digital Interpo- 
lation techniques to remotely sensed digital Images (1-3). The consensus seems 
to be that resampling using the cubic convolution method leads to the best and 
most economical results from both a subjective, statistical, and frequency con- 
tent standpoint. This attitude persists until artifacts such as ranging or other 
unpleasant visual effects are observed. The tendency at this point is to re- 
evaluate the resampling algorithm and possibly derive a new one. 

Many artifacts, which become pronounced after the application of resampling 
techniques, are present In the original data and are Independent of tl\e Interpo- 
lation method used. A resampling technique based on band limited filter analysis 
Is described in the next section, and results are compared to those obtained 
using cubic convolution on a test set of analytic functions. 

In Section III, application of this interpolation method is made to a Landsat -1 
scene which was used as a test case for the Digital Image Rectification System 
(DIRS) (4). The scene, Rockwood, Illinois, was widely discussed at the Image 
100 User's Conference held at Goddard Space Flight Center during November 
1975. Regardless of the resampling techniques used, artifacts which resemble 
ringing are present in the final image. These features are shown to be present 
in the original data and are not introduced by the x’esampling method. 


II. AN INTERPOLATION ALGORITHM 
A. General Comments 

The general criterion for interpolators or resamplers seems to be related to 
whether or not the resampled values look right and how these new values clas- 
sify. If values at edges or other locations within the image look strange, the 
tendency appears to be to change the constants and try again. Sometimes, a 
different interpolation approach may be formulated. 

On occasion one may consider the sin (x)/x as an interpolater, with the usual 
complaints about the number of side lobes required for N% accuracy and the 
possibility of introducing spurious features due to ringing. Ringing is caused 


1 


by the occurrence of sharp edges or discontinuities In the data to bo resampled 
and resembles the Interpolating function with respect to the spacing of ringing. 

It is rare tiuit one will ever encounter this situation in nature especially when 
the smoothing effects of Instrument electronics are considered. 

The process of Interpolation is fundamentally a low pass filtering operation and 
should be treated as such, especially when dealing with band limited signals. 
When a series of weights Is convolved with the picture data, the frequency 
spectrum of the lines is modified by the spectrum of the weights via multiplica- 
tion. Thus, the ringing effects of a sin (.x)/x, the Gibbs phenomenon, can be 
effectively reduced by apodlzlng processes (5). Apocllzing consists of multiply- 
ing the filter weights by a function which is unity and has zero slope at the center 
(maximum weight) and which appi’oaches zero wltli zero slope at the final weight. 

The following criterion will be used as a basis for choosing an Interpol.atlon 
function (or set of coefficients): 

Does the function give correct numeric values for any analj^tic input function at 
values between equally spaced points assuming the minimum Nyqulst sampling 
rate? (Note that oversampling will always aid us.) The rationale for this defini- 
tion is tliat over a small number of Image points along or perpendicular to the 
flight line, some "best fit" function can be defined such as an N*'* oixler poly- 
nomial. Admittedly, original data points may be altered slightly, but one can 
then generate another polynomial based on the differences at the data points, 
subtract it from the initial polynomial and still remain with an analytic ex- 
pression. 

In view of the previous discussion, the function 
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sin (^n k) 
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11 " 0, *1, +2, ... iN 


( 1 ) 


is submitted as an interpolation constant generator where 

Ic = the number of interpolated intervals between two original sample 
points. 

N = j ' k + 1 where j defines the number of zero crossings of the sin (x)/x 
side lobes along a positive or negative axis. 

This expression is the apodized Fourier Transform of an idealized low pass 
filter (5). 
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Column 1 of Table 1 tabulates this function for It = 4 and J ~ 3. If the values 
are applied to a "picket fence", that Is, equally spaced pulses of unit amplitude, 
a certain amount of ringing will occur as can be seen by summing the appropri- 
ate coefficients used In the convolution. This effect Is due to using a finite 
number of interpolation constants to approximate the picket fence at Intermedi- 
ate sample values. Normaliz,lng each Interpolation constant to the sum of the 
values Involved In the convolutloni 


j-1 

E Vk + i) 

11== -j 


1=1,2, . . . , k - 1 


( 2 ) 


effectively eliminates this effect and values are tabulated In Column 2, Table 1. 

Figure 1 Illustrates the application of the latter coefficients to a discretely 
sampled input for a midpoint interpolation. When the coefficients are slid to the 
right an amount equal to one quarter of original sampling unit, a new value may 
be calculated. In all cases ojily 6 data points and the proper choice of G of the 
interpolation coefficients are required. Note that the resampled value is located 
at the peak of Interpolation coefficients. 

Table 2 applies the Interpolation constants of Table 1, Column 2, to the discrete 
polynomials (F (n) ) designated in the headings. The values for n, the theoretical 
value of F (n), and the interpolated values are shown for each polynomial. The 
Interpolated values are calculaccd from F(n) with n being an Integer. It should 
be noted that for values midway between the original samples, the interpolated 
values are nearly exact. This leads one to consider a recursive interpolation 
technique in which values are calculated for the midpoints, and this new set of 
data is then used as a basis for interpolating new midpoint values which corre- 
sponds to computing values at one fourth of the original spacing. This process 
is illustrated In the third major heading of Table 2. F'(n) and F"(n) refer to the 
first and second passes of the interpolation process. 

Table 3 Illustrates the midpoint evaluation on the cubic polynomial as designated 
by cubic convolution metliods. As a comparison, midpoint convolution for the 
case where j = 2 and k - 2 utilizing equation 1 is shown. In both these cases, 
four data points are used thu.s agreement with the correct midpoint value is 
not as good as when six data points are used as in Table 2. Values for integer 
values of n are shown once. 

The Fourier Transforms of the function used (k = 4, j = 3) are shown in Fig- 
ure 2 for the unapodized, apodlzed, and normalized cases. The deviation from 
the ideal low pass box car is due to the small number of lobes used. Figure 3 
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shows the transCornis of upodlzcd sine functions for 4, 8 and 16 side lobes on 
either side of the maximum and the improveniont In approaching the box car Is 
noticeable- This Is nothing more than a restatement of the fact that interpolated 
values improve as more data Is used. 


III. ROCK WOOD ILLINOIS SCENE 

Figure 4 Is a DIRS (4) processed scene of the Rockwood Illinois area and corre- 
sponds to 10 meter resolution. This scene was the basis of much discussion 
during the Image 100 User's Conference held at Goddard Space Flight Center in 
the latter part of 1975. The item discussed was the striping effects observed 
in tlio lower left and center left regions of the picture. These effects were at- 
tributed to the resampling procedure which In this case was based on cubic 
convolution. 

This is not the case. The "striping" is a manifestation of the original data and 
Is due to tile way the Landsat Multlspectral Scanner (MSS) samples nai'row fea- 
tures sufficiently different from the adjacent background. As an example*, 
consider a cement road 40 meters wide surrounded by water on both sides in a 
direction which parallels the flight path. The Instantaneous field of view of the 
MSS is about 80 motors square. Due to scanner sampling timing (4), a 22 meter 
offset occurs between detectors 1 and 6 of a particular band, and earth rotation 
contributes about a .50 meter offset in opposite direction between detector 6 of 
swath n, and detector 1 of swath n + 1. Given this situation, it is easy to visu- 
alize two adjacent sp.mples from detector 1, swath n acquiring half the road for 
a net Intensity of I, and anottier sample from swath n, detector 6 acquiring the 
full road for an intensity of 21. Thus the modulation of this feature can be 50%, 
and any interpolation of resampling algorithm should reflect this. 

Figure 5 is the original data from the lower half of the scene previously de- 
scribed and consists of 150 pixels and 150 lines. Each pixel and line has been 
replicated 8 times In order to simulate 10 meter resolution. Additionally each 
line has been shifted to the nearest fraction of a pixel to allow for the MSS 
sampling skew. In the regions where Figure 3 shows striping, one observes 
that the replicated data also shows a striping tendency. The I’esampling proce- 
dure simply smoothes the transitions in these areas and appears as low frequency 
striping In Figures 4 and G. 


*'fhis analony is valid if llie sciimiiiig ruiiolloii is rcclang iliir and sees only one “pixel”. In practice die 
scanning funclion is a two dimensional smear of the order of 4 pixels by 3 linos. A deconvolution method 
based on a priori knowledge of the scanner is under investigation by Bcndlx Coiporation. 
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Flguro G was generated from the original data for Figure 4 using the recursive 
interpolation discussed in the previous section. Initially, 3 Itorations for eacli 
lino wore made (a one eightii pixel interval resampling). The linos were shifted 
to allow for scanner skew, and then three more recursive interpolations wore 
computed In a "vertical" direction co !,’responding to a one eighth line resampling. 
As can bo seen, the regions of interi st, except for scale, are identical. 

IV. CONCLUSIONS 

A band limited, convolutional, interpolation technique using few coefficients was 
developed in Section I. The accuracies obtained wlion using analytic Inputs is 
very high for new samples calculated at midpoint spacing. The Interpolated 
values in the cases illustrated are closer to the actual calculated values than 
those obtained using cubic convolution methods. 

An additional advantage of this method Is the ease with wliich the rate of cutoff 
at the cutoff frequency may bo varied by varying the number of side lobes. This 
ean be significant if there are major liigh frequency contributions in the data to 
be Interpolated. The three side Iqbe coefficients perform adequately when the 
highest frequency is sampled three or more times. When the sampling rate ap- 
proaches the Nyquist rate (2 samples per highest frequency) more side lobes are 
required to improve the filter cutoff (see Figures 2 and 3). 

Data artifacts which resemble ringing were shown to be present in the data and 
not caused by interpolation for a particular Landsat scene discussed in Section n. 
Recursive convolution, which is a natural offshoot of the bandlimited interpola- 
tion, was applied to tlie data for this scene and results were "eyeball" identical 
to previous processing. 


V. POSSIBLE FUTURE APPLICATION 

One of the initial goals of this study was to resample an image using an Inter- 
polative nearest neighbor technique. Table lookup procedures would be utilized. 
Thougli not accomplished, this could be implemented by deriving Interpolation 
coefficients for each . 1 or . 01 of a pixel, computing the proper nearest neighbor 
location and applying the looked up coefficients to the data. 

No attempts to classify results obtained with the interpolation scheme discussed 
have been ti’ied, however if "mixels" are avoided, results should be favorable. 
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Figure 2. Transfer Function of Interpolation Coefficients 




Figure 3. Approximation to an Ideal "Box Tar" as a Function of the Number of Side Lobes 





Figui’c 5. Original Rockwood, Illinois Data Replicated 8 Times 
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Figure (5. Rockwowl, Illinois Data Interpolatocl to 10 Meter Resolution by 

Recursive (’onvolution 
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Table 1 

Interpolation Function Tabulation 


n 

1 

2 

0 

1.00000 

1.00000 

± 1 

.88970 

.8925 

± 2 

.60680 

.6110 

± 3 

.26900 

.2698 

± 4 

.00000 

.00000 

± 5 

-.13070 

-.1311 

± 6 

-.13140 

-.1323 

± 7 

1 

-.06485 ' 

-.06505 

± 8 

.00000 

.00000 

i 9 

.02712 

.02721 

± 10 

.02122 

.02137 

± 11 

.00660 

.00662 

± 12 

.00000 

.00000 


Col. 1 = [l - __ 

•- K n/4 ‘ 


Col. 2 " + Is ^ ^ ^-3 ^ ^-7 ^ n) 

for n = 1, 5, 9, -3, -7, “11 etc. 
Denominator = C j 






Sample Computations Using Table 1 Coefficients 
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Note: If coefficients are applied to 1, -1, 1, -1, — , 1, -1 in recursive manner sin{x) will be calculated to 
bstter*than 3% at values of t/2" 




Table 3 


Comparison of Interpolation 
Using Cubic Convolution and Interpolation Function 

F (n) = n“ 8000 


n 

F(n) 

Cubic Conv. 
F(n) 


18 

-2168 



19 

-1141 



20 

0 



20.5 

615.1 

600 

609 

21 

1261 



22 

2648 



23 

4167 



23.5 

4977.9 

4960 

4969 

24 

5824 



25 

7625 
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